In [ ]:
from __future__ import division, print_function
import os
import sys
from six.moves import cPickle as pickle
# Third-party
import astropy.coordinates as coord
import astropy.units as u
uno = u.dimensionless_unscaled
import matplotlib as mpl
import matplotlib.pyplot as pl
%matplotlib inline
pl.style.use('apw-notebook')
import numpy as np
from scipy.signal import argrelmax, argrelmin
# Custom
importgala.coordinates as gc
importgala.dynamics as gd
fromgala.dynamics.util import estimate_dt_nsteps
importgala.integrate as gi
importgala.potential as gp
fromgala.units import galactic
import superfreq
from ophiuchus import RESULTSPATH
import ophiuchus.potential as op
from ophiuchus.data import OphiuchusData
from ophiuchus.plot import plot_data_orbit
from ophiuchus.experiments import LyapunovGrid
plotpath = "/Users/adrian/projects/ophiuchus-paper/figures/"
In [ ]:
ophdata = OphiuchusData()
In [ ]:
all_names = ["static_mw"] + ["barred_mw_{}".format(i) for i in range(1,10)]
short_names = ["static"] + ["bar{}".format(i) for i in range(1,10)]
name_map = dict(zip(all_names, short_names))
In [ ]:
# name = 'barred_mw_1'
# gr = LyapunovGrid.from_config(os.path.join(RESULTSPATH,name,"lyapunov"),
# os.path.join(RESULTSPATH,"global_lyapunov.cfg"),
# potential_name=name)
# d = gr.read_cache()
# 1/d['mle_end']/1000.
In [ ]:
name = 'static_mw'
# bins = np.linspace(0.3,1.1,12)*u.Gyr
gr = LyapunovGrid.from_config(os.path.join(RESULTSPATH,name,"lyapunov"),
os.path.join(RESULTSPATH,"global_lyapunov.cfg"),
potential_name=name)
d = gr.read_cache()
ftmle = (d['mle_end']*1/u.Myr)
lyap_time = (1/ftmle).to(u.Myr)
pl.hist(lyap_time, bins=16)
In [ ]:
fig,axes = pl.subplots(3,3,figsize=(7,7), sharex=True, sharey=True)
bins = np.linspace(0.3,1.1,12)*u.Gyr
for i,name in enumerate(all_names[1:]):
gr = LyapunovGrid.from_config(os.path.join(RESULTSPATH,name,"lyapunov"),
os.path.join(RESULTSPATH,"global_lyapunov.cfg"),
potential_name=name)
d = gr.read_cache()
ftmle = (d['mle_avg']*1/u.Myr)
lyap_time = (1/ftmle).to(u.Myr)
axes.flat[i].hist(lyap_time, bins=bins.to(u.Myr))
# axes.flat[i].set_title(name_map[name], fontsize=18)
axes.flat[i].text(1100, 45, name_map[name], fontsize=18, ha='right')
# axes.flat[i].axvline(lyap_time[0].value, color='#2b8cbe', linestyle='dashed', lw=2., ymax=37/50*1)
# if i > 5:
# axes.flat[i].set_xlabel(r"$t_\lambda$ [Myr]", fontsize=18)
axes[0,0].set_xlim(300,1200)
axes[0,0].xaxis.set_ticks([400,750,1100])
axes[0,0].yaxis.set_ticks([0,20,40,60])
axes[1,0].set_ylabel('$N$')
axes[2,1].set_xlabel(r"$t_\lambda$ [Myr]", fontsize=18)
axes[0,1].set_title(r"Larger pattern speed $\longrightarrow$", fontsize=18, y=1.04)
axes[1,2].set_ylabel(r"$\longleftarrow$ Larger bar angle", fontsize=18, labelpad=10)
axes[1,2].yaxis.set_label_position("right")
# fig.savefig(os.path.join(plotpath, "lyapunov-hist.png"), dpi=300)
# fig.savefig(os.path.join(plotpath, "lyapunov-hist.pdf"))
In [ ]:
all_lyap_times = np.array([])
all_periods = np.array([])
all_Omega = np.array([])
for i,name in enumerate(all_names[1:]):
print(name)
# integrate and estimate periods
w0 = np.load(os.path.join(RESULTSPATH,name,"orbitfit","w0.npy"))[:127]
pot = op.load_potential(name)
orbits = pot.integrate_orbit(w0.T, dt=-1, nsteps=1000, Integrator=gi.DOPRI853Integrator)
periods = np.abs([orbits[:,i].estimate_period().value for i in range(127)])*u.Myr
# read lyapunov times
gr = LyapunovGrid.from_config(os.path.join(RESULTSPATH,name,"lyapunov"),
os.path.join(RESULTSPATH,"global_lyapunov.cfg"),
potential_name=name)
d = gr.read_cache()
ftmle = (d['mle_avg']*1/u.Myr)
lyap_times = (1/ftmle).to(u.Myr)
# color by pattern speed
Omega = np.zeros_like(periods.value) + pot.parameters['bar']['Omega']
all_lyap_times = np.concatenate((all_lyap_times, lyap_times.value))
all_periods = np.concatenate((all_periods, periods.value))
all_Omega = np.concatenate((all_Omega, Omega))
In [ ]:
pl.hist(all_lyap_times, bins=np.linspace(400,1500,32))
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(7,7))
c = ax.scatter(all_lyap_times, all_periods, c=all_Omega)
fig.colorbar(c)
pl.xlim(400,1000)
# axes.flat[i].text(1100, 45, name_map[name], fontsize=18, ha='right')
# axes.flat[i].axvline(lyap_time[0].value, color='#2b8cbe', linestyle='dashed', lw=2., ymax=37/50*1)
# if i > 5:
# axes.flat[i].set_xlabel(r"$t_\lambda$ [Myr]", fontsize=18)
# axes[0,0].set_xlim(300,1200)
# axes[0,0].xaxis.set_ticks([400,750,1100])
# axes[0,0].yaxis.set_ticks([0,20,40,60])
# axes[1,0].set_ylabel('$N$')
# axes[2,1].set_xlabel(r"$t_\lambda$ [Myr]", fontsize=18)
# axes[0,1].set_title(r"Larger pattern speed $\longrightarrow$", fontsize=18, y=1.04)
# axes[1,2].set_ylabel(r"$\longleftarrow$ Larger bar angle", fontsize=18, labelpad=10)
# axes[1,2].yaxis.set_label_position("right")
# fig.savefig(os.path.join(plotpath, "lyapunov-hist.png"), dpi=300)
# fig.savefig(os.path.join(plotpath, "lyapunov-hist.pdf"))
In [ ]:
In [ ]:
import time
In [ ]:
nperiods = 1024
nsteps_per_period = 1024
fac = nsteps_per_period/256.
t0 = time.time()
lyap = gd.fast_lyapunov_max(gr.w0[0], op.load_potential('barred_mw_1'),
dt=d['dt'][0]/fac, nsteps=nperiods*nsteps_per_period,
return_orbit=False)
print(time.time()-t0)
In [ ]:
pl.loglog(np.arange(lyap.shape[0])/fac, 1/lyap/1000., marker=None)
In [ ]:
pl.loglog(np.arange(lyap.shape[0])/fac, 1/lyap/1000., marker=None)
In [ ]:
name = 'barred_mw_fixed'
with open(os.path.join(RESULTSPATH,name,"lyapunov","orbit.pickle"),'rb') as f:
orbit = pickle.load(f)
In [ ]:
potential = op.load_potential(name)
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(4,4))
fig = potential.plot_contours(grid=(np.linspace(-15,15,32),np.linspace(-15,15,32),0.), ax=ax)
In [ ]:
fig = orbit[:,0].plot(linestyle='none', marker=',', alpha=0.1)
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(5,5))
with open(os.path.join(RESULTSPATH,name,"lyapunov","lyap.pickle"),'rb') as f:
lyap = pickle.load(f)
lyap = np.mean(lyap, axis=1)
# ax.plot(ts[10:-10:10], 1/lyap, marker=None, color=color, label=label)
ax.plot((1/lyap).to(u.Gyr), marker=None)
ax.set_xlabel("$N$ iterations")
ax.set_ylabel(r"$t_\lambda$ [Gyr]")
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(10,lyap.size)
ax.set_ylim(1E-2,22)
fig.tight_layout()
In [ ]:
# w0 = orbit[0,0]
w0 = gd.CartesianPhaseSpacePosition(pos=[8.,0,0]*u.kpc, vel=[0,220,20.]*u.km/u.s)
f_orbit = potential.integrate_orbit(w0, dt=2., nsteps=64000, Integrator=gi.DOPRI853Integrator)
w = f_orbit.w(units=galactic)
t = f_orbit.t.value
In [ ]:
tt = t[:f_orbit.t.size//2+1]
sf = superfreq.SuperFreq(tt)
fs1 = [(w[i,:t.size//2+1,0]+1j*w[i+3,:t.size//2+1,0]) for i in range(3)]
fs2 = [(w[i,t.size//2:,0]+1j*w[i+3,t.size//2:,0]) for i in range(3)]
res1 = sf.find_fundamental_frequencies(fs1, min_freq_diff=1E-4)
res2 = sf.find_fundamental_frequencies(fs2, min_freq_diff=1E-4)
In [ ]:
res1.fund_freqs, res2.fund_freqs
In [ ]:
frac_rate = (res1.fund_freqs - res2.fund_freqs) / res1.fund_freqs / t[t.size//2]
In [ ]:
frac_diff_time = (1 / frac_rate) * u.Myr
frac_diff_time.to(u.Gyr) / 1000. # globular cluster scale
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(5,5))
colors = ['k'] + ['#B2182B']*9
# labels = ['static', 'barred'] + [None]*8
labels = [None]*10
for name,color,label in zip(all_names,colors,labels):
with open(os.path.join(RESULTSPATH,name,"lyapunov","lyap.pickle"),'rb') as f:
lyap = pickle.load(f)
lyap = np.mean(lyap, axis=1)
# ax.plot(ts[10:-10:10], 1/lyap, marker=None, color=color, label=label)
ax.plot((1/lyap).to(u.Gyr), marker=None, color=color, label=label)
ax.set_xlabel("$N$ iterations")
ax.set_ylabel(r"$t_\lambda$ [Gyr]")
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(10,lyap.size)
ax.set_ylim(1E-2,22)
ax.text(1.2E4, 10, "Static", rotation=50, fontsize=18)
ax.text(1.2E4, 2.5E-1, "Barred", fontsize=18, color=colors[-1])
# ax.legend(loc='upper left', fontsize=18)
# ax.xaxis.set_ticklabels(["", "1", "10", "100", "1000", "10000"])
# ax.yaxis.set_ticklabels(["", "0.1", "1", "10", "100", "1000"])
fig.tight_layout()
# fig.savefig(os.path.join(plotpath, "lyapunov.png"), dpi=300)
# fig.savefig(os.path.join(plotpath, "lyapunov.pdf"))
In [ ]:
ts = dict()
ws = dict()
for name in all_names:
w0 = np.load("/Users/adrian/projects/ophiuchus/output/orbitfit/{}/w0.npy".format(name))
w0 = np.median(w0, axis=0)
pot = op.load_potential(name)
dt, nsteps = estimate_dt_nsteps(w0, pot, nperiods=256, nsteps_per_period=256, dE_threshold=None)
ts[name],ws[name] = pot.integrate_orbit(w0, dt=dt, nsteps=nsteps, Integrator=gi.DOPRI853Integrator)
# testing rotating coords
xs = ws[name][:,0,:3]
fig = gd.plot_orbits(xs, marker=None)
In [ ]:
freqs = np.zeros((len(all_names),2,3))
for j, name in enumerate(all_names):
t = ts[name]
w = ws[name]
tt = t[:t.size//2+1]
sf = superfreq.SuperFreq(tt)
fs1 = [(w[:t.size//2+1,0,i]+1j*w[:t.size//2+1,0,i+3]) for i in range(3)]
fs2 = [(w[t.size//2:,0,i]+1j*w[t.size//2:,0,i+3]) for i in range(3)]
freq1,tbl1,ixes1 = sf.find_fundamental_frequencies(fs1, min_freq_diff=1E-4)
freq2,tbl2,ixes2 = sf.find_fundamental_frequencies(fs2, min_freq_diff=1E-4)
# vecs1 = superfreq.find_integer_vectors(f1, tbl1)
# vecs2 = superfreq.find_integer_vectors(f2, tbl2)
f1 = [np.sum(tbl1[tbl1['idx']==k]['A'][None] * np.exp(1j * tbl1[tbl1['idx']==k]['freq'][None] * tt[:,None]), axis=1) for k in range(3)]
f2 = [np.sum(tbl2[tbl2['idx']==k]['A'][None] * np.exp(1j * tbl2[tbl2['idx']==k]['freq'][None] * tt[:,None]), axis=1) for k in range(3)]
freqs[j,0] = freq1
freqs[j,1] = freq2
In [ ]:
fig,axes = pl.subplots(2,1,figsize=(6,8),sharex=True,sharey=True)
axes[0].plot(tt, fs1[1].real, marker=None)
axes[1].plot(tt, f1[1].real, marker=None)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
reg_w0 = np.array([25.,0,0,0.,0.01,0.2]) # regular
oph_w0 = w0[0] # ophiuchus
In [ ]:
# reg_t,reg_w = pot.integrate_orbit(reg_w0, dt=1., nsteps=250000, Integrator=gi.DOPRI853Integrator)
oph_t,oph_w = pot.integrate_orbit(oph_w0, dt=0.4, nsteps=250000, Integrator=gi.DOPRI853Integrator)
In [ ]:
fig = gd.plot_orbits(reg_w, marker=None)
fig = gd.plot_orbits(oph_w, marker=None, axes=fig.axes)
In [ ]:
# Lz = np.cross(w[:,0,:3], w[:,0,3:])[:,2]
# EJ = pot.total_energy(w[:,0,:3], w[:,0,3:]) - Omega * Lz
In [ ]:
# dE = np.abs((EJ[1:]-EJ[0])/EJ[0])
# pl.semilogy(dE, marker=None)
In [ ]:
for t,w in zip([reg_t,oph_t], [reg_w[:,0], oph_w[:,0]]):
sf = superfreq.SuperFreq(t[:t.size//2])
fs1 = [(w[:t.size//2+1,i]+1j*w[:t.size//2+1,i+3]) for i in range(3)]
fs2 = [(w[t.size//2:,i]+1j*w[t.size//2:,i+3]) for i in range(3)]
f1,_,_ = sf.find_fundamental_frequencies(fs1)
f2,_,_ = sf.find_fundamental_frequencies(fs2, min_freq_diff=1E-4)
print(f1)
print(f2)
print(t.max() / (2*np.pi/f1))
print("--")
In [ ]:
logpot = gp.LogarithmicPotential(v_c=1., r_h=0.1, q1=1., q2=0.8, q3=0.8, units=galactic)
Omega_b = np.array([0.,0.,0.1])
def func(t, w):
q = w.T[:3].T
p = w.T[3:].T
dq = p - np.cross(Omega_b[None], q)
gradPh = logpot.gradient(q)
dp = -gradPh - np.cross(Omega_b[None], p)
return np.hstack((dq, dp))
In [ ]:
integrator = gi.DOPRI853Integrator(func)
In [ ]:
t,w = integrator.run(np.array([1.,0,0,0,-0.8,0]), dt=0.1, nsteps=10000)
In [ ]:
pl.plot(w[:,0,0], w[:,0,1])
In [ ]: